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Abstract. We study the pumping of magnetic flux in three-dimensional compressible magnetoconvection in the 
context of stellar dynamos. The simulation domain represents a rectangular section from the lower part of a stellar 
convection zone plus the underlying stably stratified layer, with a total depth of up to five pressure scale heights. 
Once convection has attained a statistically stationary state, a magnetic field is introduced. The magnetic field 
is subsequently modified by the convective motions, and the resulting pumping effects are isolated by calculating 
various coefficients of the expansion of the electromotive force, u x b, in terms of components of the mean magnetic 
field. The dependence of the pumping effects on rotation, latitude and other parameters is studied. First numerical 
evidence is found for the existence of pumping effects in the horizontal directions, unless the rotation axis coincides 
with the vertical axis, as is the case on the poles. Evidence is found that the pumping effects act differently on 
different components of the mean magnetic field. Latitudinal pumping is mainly equatorward for toroidal field, 
and can be poleward for poloidal field. Longitudinal pumping is mainly retrograde for the radial field but prograde 
for the latitudinal field. The pumping effect in the vertical direction is found to be dominated by the diamagnetic 
effect, equivalent to a predominating downward advection with a maximum speed in the turbulent case of about 
10% of the rms convective velocity. Where possible, an attempt is made to identify the physical origin of the 
effect. Finally, some consequences of the results for stellar dynamos are discussed. 
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1. Introduction 

Magnetic pumping refers to transport of magnetic fields 
in convective layers which does not result from bulk mo- 
tion. The existence of magnetic pumping has been known 
for several decades now. Drobyshevski & Yuferev (1974) 
demonstrated in a numerical investigation that magnetic 
flux in a convective layer is expelled downward if the down- 
flows form connected lanes that separate isolated upflow- 
ing regions, an effect known since as topological pump- 
ing (cf. also Moffatt [1978| , p. 70; Petrovay |l99l| ) . A dif- 
ferent but related effect is turbulent pumping. For the 



(Kichatinov & Pipin 1993). Numerical simulations of 



stratified convection showed that dynamo-generated mag- 



netic fields are pum ped d ownward (Nord lund e t al. |1992 
Brandenburg et al. |1996[) Tobias et al. ( |l99S| , [200 1|) and 



two-dimensional case, Drobyshevski (1977) showed that 
a poloidal magnetic field is transported against the gra- 
dient of the density. Analytical investigations based on 
the first order smoothing approximation (FOSA) suggest 
that turbulent pumping also occurs in the presence of a 
gradient in the intensity of the turbulence, due t o the 
diamagnetic effec t (Ra dlcr 1968 , Roberts & Soward 1975J, 



Dorch & Nordlund ( 2001 ) performed detailed numerical 
experiments of magnetic pumping in box simulations of 
magnetoconvection. They did not, however, measure the 
pumping velocity which plays a role in mean-field electro- 
dynamics. 

Although in most investigations the possible signifi- 
cance of magnetic pumping for stellar dynamo action is 
stressed, it was in fact rarely incorporated in models. This 
is true partly because early mean-field solar dynamo mod- 
els did not appear to 'need' any pumping effect since 
they were able to reproduce the large-scale solar mag- 
netic field rather well using the simplest assumption that 
there are just alpha effect and diffe rential rotation (for 
the alpha effect see, e.g., Parker 195E, or Krause & Radlcr 



Krause & Radler |1980| , Kichatinov & Riidiger |1992[ ) , or as 
a result of rotationally induced anisotropies (Kichatinov 
1991), or through magnetic buoyancy of small-scale fields 



1980 ). Hence the pumping effect, along with many other 



dynamo coefficients, was thought to be an unnecessary 
complication. First indications for the usefulness of tur- 
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bulent pumping came from an af2-dynamo simulation of 
Brandenburg et al. ( I992j ), where downward pumping of 
the mean magnetic field was found to favour an equator- 
ward migration of magnetic field in the presence of mainly 
latitudinal differential rotation. 

A major problem in solar dynamo theory was the 
storage of magnetic flux for suffi ciently long times (see, 
e.g., Moreno-Insertis et al. 1992 ). It was believed that 



magnetic buoyancy would remove magnetic field quickly 
from the convection zone. However, numerical simulations 
of dynamo action in turbulent stratified convection indi- 
cated quite clearly that the effects of magnetic buoyancy 
can be offs et by the downward p umpi ng (Brandenburg & 
Tuomincn 1991, Nordlund et al. 1992, Brandenburg et al. 



I99E, Sect. 4.4). Thus the vertical pumping effect may be 



important because sufficient amplification of the toroidal 
magnetic fields generated by the solar dynamo can be 
achieved only in or near the stably stratified overshoot 
layer below the convection zone. 

Apart from the vertical direction, pumping effects 
can also exist in other directions. So far though, non- 
radial pumping effects have not received much attention, 
although they have interesting consequences for stellar 
dynamos, as is discussed further below. Their existence 
has been predicted on the basis of symmetry considera- 
tions within the framework of mean-field electrodynamics 



(Krause & Radler 1980), and some analytical results based 
on FOSA are presented in Kichatinov ( I991| ) and Rudiger 
& Kichatinov ( 1993 ). No numerical studies of non-radial 
pumping effects have previously been carried out, as far as 
we know. The present investigation provides the first nu- 
merical evidence for their existence, as well as quantitative 
results for the corresponding transport coefficients. 

The study of pumping effects is part of an effort to 
calculate dynamo coefficients through numerical simula- 
tions. This work is motivated by the idea that progress 
in understanding stellar dynamo action can be achieved 
only through a combination of approaches which include 
numerical simulation as well as mean-field modelling. On 
the one hand, fully global MHD calculations of stellar 
dynamos have not been successful so far because com- 
puter capacities are insufficient to capture the parameter 
regime of stellar convection zones. Local MHD simulations 
of small subvolumes offer the best possibility to approach 
solar parameters. On the other hand, observations of the 
large-scale solar magnetic field (solar cycle, butterfly di- 
agram, Hale's polarity laws) suggest that some form of 
mean field description is possible. Local numerical inves- 
tigations of aspects of stellar dynamos may provide useful 
input information for the construction of mean-field mod- 
els. Also, the study of dynamo coefficients can be a pow- 
erful diagnostic tool to isolate from numerical simulations 
the various effects that play a role in stellar dynamos and 
to quantify them. In order to understand the role of dy- 



namo coefficients, consider the averaged MHD induction 
equation: 



OB 

— = curl (UxB + S -77 curl B) , 



(1) 



Here B is the mean magnetic field, U is the mean flow, rj 
is the magnetic diffusivity, and S = u x b is the electro- 
motive force (EMF) which depends on u = U — U, the 
fluctuating component of the flow, and b = B — B, the 
fluctuating component of the magnetic field. Mean quan- 
tities are either ensemble averages or azimuthal averages. 
In the present case of a box with periodic boundary con- 
ditions in the horizontal directions (see below), horizontal 
averages can be used instead, and in practice, there is ad- 
ditional temporal averaging. The dynamo coefficients ap- 
pear in an expansion of the electromotive force in terms of 
spatial derivatives of the mean magnetic field. In Cartesian 
coordinates, and assuming the summation convention, this 
may be written as 



Si = otijBj 



PijkV jBk 



(2) 



The coefficients in the expansion, of which only the first 
two are shown, are pseudo tensors of increasing rank. In 
what follows, we consider only the first term, i.e. we as- 
sume that the mean magnetic field in the simulation do- 
main is approximately independent of spatial coordinates. 
Commonly one separates the a-tensor into symmetric and 
antisymmetric parts, i.e. 



(3) 



The latter gives rise to an advection term, a A ■ B = 7 x B, 
where 



o 

7=- I a 



A 

zx 



(4) 



This term describes the pumping of all components of the 
magnetic field with equal speed, and we shall refer to it as 
the 'general pumping effect'. Note that the pumping effect 
as it is used here refers to the zeroth-order terms in the 
expansion of the EMF only. Turbulent diffusion ((3) and 
higher order contributions can also give rise to transport 
phenomena of the mean magnetic field, but they cannot be 
expressed by advection terms and are to be treated sepa- 
rately. One may also ask whether there are circumstances 
under which different components of the magnetic field 
have different advection speeds and directions. A neces- 
sary condition that must be fulfilled when addressing this 
question is that the magnetic field can be decomposed 
in divergence-free components. Only then can one write 
down separate equations for the magnetic flux of each field 
component through a surface moving with a velocity w, 
say, thereby allowing one to distinguish advection from 
other processes. This holds, for example, for a decomposi- 
tion using toroidal and poloidal components (Kichatinov 
1991 ). But unlike Alfven's theorem in ideal MHD (where 
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w = U for all divergence- free components), it is gener- 
ally impossible in mean-field electrodynamics to define a 
moving surface in which flux is conserved. That is due to 
the presence of the alpha effect and turbulent diffusion. 
Essentially, what can be achieved is to define a velocity 
that optimizes flux conservation for a given divergence- 
free component. For example, in the case of B t0 r, one sets 
Wtor = f7+7' tor ' , where 7( tor ) denotes the pumping veloc- 
ity of the toroidal field (even though it is itself poloidal). 

As a result of the averaging over the horizontal co- 
ordinates, our mean magnetic field is of the form B = 
[B x (z), B y (z), B z ], so that all three Cartesian compo- 
nents are divergence free. If one interprets this field as 
the mean field of an axisymmetric spherical dynamo, then 
B x e x + B z e z corresponds to the poloidal component and 
B y e y to the toroidal. The meaning of the tensor compo- 
nents a.ij may be grasped most easily by writing 



B 



j {x) x B x +7 (a) x By 



7 {z) xB z , 



(5) 



where B x = B x e x , etc. Here, a D denotes the diagonal 
part of a, which is responsible for maintaining a large- 
scale magnetic field in solar-type stars. The off-diagonal 
terms of the a tensor represent various contributions to 
the pumping effect. As shown by Eq. (|J), they are equiv- 
alent to an effective advection of the mean magnetic field, 
where 



7 



7 



7- 



7 












(G) 



are the pumping velocities for the respective Cartesian 
components of the mean field. This representation shows 
that the pumping velocity can, in principle, be different for 
each component, as was pointed out by Kichatinov (1991 ) , 
and that the differences are equal to the off-diagonal, sym- 
metric components of the a-tensor. 

In general, the nine components are independent. 
At the poles however, the rotation axis coincides with 
the direction of gravity, resulting in isotropy in the hori- 
zontal plane. Since the only remaining preferential direc- 
tion is the vertical one, it follows from elementary sym- 
metry considerations that the pumping velocity cannot 
have horizontal components, i.e. a xz — a zx — a yz — 
&zy = 0, and that the pumping effect acts in identical 



ways on B x and B y: i.e. a x 



-a yx . By the same to- 



ken a xx — a yy . Hence at the poles the EMF is given by 
£ = a xx Bn + a zz B z + 7 x B, where B H = [B x , B v , 0] is 
the horizontal mean field and 7 = —a xy e z is the pump- 
ing velocity. At the equator, similar arguments lead to the 
conclusion that the pumping velocity in the x-direction 
vanishes, i.e. a yz — a zy — 0, and that a D = 0. These 
same components are antisymmetric with respect to the 
equatorial plane, while the other ones are symmetric. 




Fig. 1. Geometry of the simulation domain. 

The remainder of the paper is structured as follows. 
After introducing the model (Sect. 2), we explain the 
method by which the dynamo coefficients are calculated 
(Sect. 3). In Sect. 4 we present results on the pump- 
ing effects, their dependence on various parameters, and 
we make a qualitative comparison with analytical results, 
where they are available. Briefly, some results on the di- 
agonal alpha coefficients will be mentioned. Finally, the 
results and their possible implications for stellar dynamos 
are discussed (Sect. 5). 

2. Model 

For a detailed description of the model and for the full 
equations and boundary conditions, the reader is referred 
to Ossendrijver et al. ( |2001 ; hereafter cited as Paper I). 
The simulation domain consists of a rectangular box which 
is defined on a Cartesian grid, with x and y denoting the 
two horizontal coordinates (corresponding to latitude and 
longitude in spherical coordinates), and z denoting depth, 
respectively (Fig. Q). Gravity is directed in the positive in- 
direction. The box rotates about an axis in the xz-plane, 
Q, — (e x sin 8 + e z cos 9) , where 9 is the angle between 
the rotation vector and the z-axis, or the angular distance 
to the south pole. 

From top to bottom, the simulation domain consists of 
a thin cooling layer, a convectively unstable layer, and a 
stably stratified layer with overshooting convection. The 
consecutive boundaries of the layers are located at depths 
Z\, Z2, Z3 and 2:4, where zi = and 2:3 = d in all cases. The 
horizontal extent of the box, L x L, is chosen such that 
there are always at least a few convective cells within the 
box. The governing equations are those describing mag- 
netic induction, mass continuity, and the balance of mo- 
mentum and energy. Solutions are governed by the follow- 
ing independent dimcnsionlcss parameters. The Prandtl 
number and the magnetic Prandtl number are defined as 



Pm = vjf]. 



(7) 



where xo = ^H&APn denotes a reference value of the ra- 
diative diffusivity for region 2. In order to avoid possible 
confusion with the pumping effect we write for the ratio 
of specific heats 7 a d = 5/3. The parameter £0 determines 
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the internal energy at the top, eo, and hence the pressure 
scale height there, i.e. 

£o = (7ad - l)e /gd. (8) 

The initial thermal structure of the three regions is char- 
acterized by polytropic indices m t = (1 — Vi)/Vj, where 
Vj = (d\nT/d\np)i denotes the radiative temperature 
gradient (i — 1,2,3). The adiabatic temperature gradient 
is V a( j = (7ad — l)/7ad- A measure for instability is pro- 
vided by the (initial) superadiabaticity, Si — Vj — V a d, 
which is positive in an unstably stratified medium. The 
reference Rayleigh number is defined as 

„ d qSo , . 

Ra = , 9 

where H p h = £ d + 0.5d/(m 2 + 1) is the pressure scale 
height in the middle of the unstable layer, and the pa- 
rameters S2 and Hph refer to the initial stratification. The 
Taylor number, 

Ta = {2nd 2 /vf , (10) 

measures the importance of rotation relative to viscous 
dissipation. Finally, 00 represents the rate at which inter- 
nal energy adjusts to the specified reference value eo- The 
remaining parameters are secondary. The Coriolis num- 
ber, or inverse Rossby number, measures the importance 
of the Coriolis force: 

Co = 2fir. (11) 

Here r = d/it rms is an estimate of the correlation time, 
and u rms is the rms velocity defined in a suitable way 
(e.g., by averaging over time and over the bulk of the 
convection zone in the box). In the initial state, the ra- 
diative energy flux, F ra d,2 = —ude/dz, is assumed to be 
constant throughout the domain. This determines the ra- 
diative conductivities in the three regions according to 
K1/K2 — (mi + l)/(m2 + 1). In order to obtain a con- 
tinuous function, k is smoothed across thin intermediate 
layers between the three regions. The initial stratification, 
obtained from a simplified mixing length model, is calcu- 
lated with the constraint p(zs) = po- Kinetic and magnetic 
Reynolds numbers are defined as 

Re = u rms d/ v , Rm = u lms d/ rj. (12) 

Dimensionless units are introduced by setting d = po = 
g = po = 1, where po is the vacuum permeability. 

3. Method 

Initially, a simulation is performed without a magnetic 
field until convection has attained a statistically stationary 
state. Then at most three numerical experiments are done 
with imposed weak homogeneous magnetic fields in or- 
thogonal directions. For each run, the EMF and the mean 
magnetic field are calculated, and by combining the re- 
sults, the nine coefficients are determined through the 



Table 1. Parameters. In all runs, Z\ — —0.07, Z2 = 0, 
Z3 = 1, Z4, = 1.93, L = 4, £0 = 0.1, mi = 0, 1TI2 — 1 (i.e., 
82 = 0.1), m 3 = 3, B Q = 0.001, Pr = 0.4 and Pm = 1 (i.e., 
Rm = Re). From these parameters, and Ra, v is derived 
using Eq. (^). The quantities u rms , Ma, Re and Co are 
averages over time and over the unstable layer; A lnp cz is 
the number of pressure scale heights accross the unstable 
layer. The unit of length is d, that of u rms is yfdg, and 
that of v is yj d 3 g. 
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Ta 
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10 2 
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v ■ 10 3 


0.34 


0.34 


0.17 


3.4 


3.4 


Wrms 


0.090 


0.086 


0.089 


0.11 


0.10 


Co 


2.4 


0.75 


1.2 


0.3 


2 


Ma 


0.09 


0.09 


0.11 


0.12 


0.12 


Re 


260 


260 


524 


32 


32 


Alnp cz 


3.7 


3.7 


3.8 


2.8 


2.8 



relation £ = a ■ B. A necessary condition for the results 
to be valid is that the statistical properties of the flow 
are independent of the imposed magnetic field, i.e. that 
no additional anisotropies are introduced. This is ensured 
because during the simulation, the field strength remains 
small compared to the equipartition field strength, so that 
the Lorentz force is weak. 

As was pointed out in Paper I, the determination of 
dynamo coefficients with sufficient accuracy requires both 
spatial averaging and temporal averaging over at least sev- 
eral turnover times. Compared to Paper I, we employ an 
improved method that gives more accurate results than 
can be obtained by performing an average over the hori- 
zontal coordinates plus a simple time average. The method 
aims to minimize contributions to the EMF (||) due to first 
and higher order derivatives of the mean magnetic field. 
Such contributions can arise because with the exception 
of B z , mean quantities, defined as horizontal averages, 
can have a dependence on z. Thus, a few turnover times 
after the introduction of a homogeneous magnetic field, 
the mean field is characterized by rather steep gradients 
at several locations. The only exception is B z , which by 
virtue of the periodic boundary conditions in the horizon- 
tal directions must remain unchanged. In the ideal situa- 
tion where the z-derivatives of the mean field vanish, the 
only contributing coefficients are otij , and they can be cal- 
culated exactly by performing a series of (at most) three 
runs using imposed magnetic fields with orthogonal ori- 
entations. Note that for 9 = 0° or 180° a pumping effect 
exists only in the z-direction; it can be measured from a 
single run with an imposed field in the x or y directions. In 
the presence of significant nonvanishing gradients on the 
other hand, one could account for the higher-order contri- 
butions by doing more experiments using imposed mag- 
netic fields that have gradients and calculating the corre- 
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Fig. 2. Averaging procedure as applied to run A45X. The thin curve is £ z at z — 0.34, averaged over the horizontal 
coordinates only, as a function of time (in units of (d/ g) 1 / 2 ). Vertical lines demark the times at which the magnetic 
field is initialised (T = 25), and gray bands denote the optimal subinterval, [ii,^], used for the temporal averaging. 
The resulting cumulative average, shown as dots connected by thick straight lines, converges after typically a few 
initialisations. For the parameters of this run, see Table FD. 
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Fig. 3. Averaging procedure as applied to run A45X. Top: 
£ z as a function of time at 9 equidistant locations between 
z = —0.07 (top curve) and z = 1.55 (bottom curve). A 
vertical shift has been applied for clarity; all curves be- 
gin at the origin. The dotted vertical lines demark the 
subinterval [t\, t 2 ]- For the third curve (z = 0.34), the av- 
erage over this interval corresponds to the last point in 
Fig. ||. Bottom: B x as a function of depth (in units of 
d) at t = 0,4,8,12,24. At t = 0, the magnetic field is 
set to B x = Bq = 0.001. The averages in the lower two 
figures are calculated from Eq. (|l3|), with T = 25 and 
N = 21. The vertical extent of the shaded areas is 2cr mcan 
(see text). For the parameters of this run, see Table [|. 



sponding dynamo coefficients. Even from one experiment 
one can calculate the additional diffusive contributions if 
the mean field is varying in space and time. This allowed 



Brandenburg & Sokoloff (2002) to estimate all components 
of the a and j3 tensors (Eq. |2j) from local simulations of 
accretion discs which feature migratory dynamo cycles. 

Instead of following this very cumbersome recipe, we 
have proceeded by periodically resetting the magnetic field 
in the entire box to the same initial homogeneous field at 
times to + iT, where i = 0, 1, • • •, N, thereby allowing the 
field to evolve for only a limited duration T. In this man- 
ner, a sequence of N time slices is obtained which, to good 
approximation, are independent realisations. By averaging 
a quantity, say A, over the horizontal coordinates and over 
the time slices according to 



A(z,t) 



1 



JV-l 



NL 2 



1 1 dx dy A(x,y, z,to +t + iT), (13) 



where < t < T, one obtains something equivalent to an 
ensemble average of a time slice. The 'error' in A(z,t) is 
estimated from the standard deviation, er, obtained from 
all the time slices according to cr mcan = a j \f~N . In order to 
obtain significant results N must exceed a minimal value 
which depends on several factors such as the magnitude 
of the coefficient that one tries to measure, the amount 
of spatial averaging, and the level of turbulence; for the 
present runs, it was found that N should be at least 5. 

The most delicate part of the method consists of iden- 
tifying an optimal subinterval within the mean time slice, 
say t\ < t < that provides the most reliable estimate 
of the dynamo coefficients. The best estimate will be de- 
fined as the time average over this subinterval. Two con- 
ditions must ideally be fulfilled. First, the EMF must have 
attained the equilibrium value given by Eq. (|2|). This con- 
dition is not trivially fulfilled, because £ = u x b = 
at times to + iT, in apparent contradiction with Eq. (^|) . 
Obviously the paradox arises because Eq. (||) is not valid 
instantaneously from the moment when a homogeneous 
magnetic field is introduced, since it takes a while for the 
magnetic fluctuations to develop. Thus t\ must be chosen 
sufficiently late. Secondly, ti must be chosen such that 
the gradients of the mean field are still numerically small 
compared to the imposed field itself divided by, say, the 
thickness of the unstable layer. If this condition is satisfied, 
then Oiij is the dominant contribution to the EMF, pro- 
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vided that does not vanish for reasons of symmetry or 
other. Both conditions can be met if the EMF evolves suffi- 
ciently more rapidly than does the mean magnetic field. In 
practice, the time scales are rather similar in the present 
simulations, and both conditions are at best marginally 
fulfilled. 

The method is illustrated in Figs. |^ and |3| for run 
A45X (see Table [l] for the parameters of this run) . Fig. ^| 
shows the z-component of the EMF at a fixed depth in 
the box, defined as an average over x and y, as a function 
of time. Although the temporal evolution varies strongly 
between the different initialisations, the additional aver- 
aging over the subinterval [ti,ta] an( i over the subsequent 
realisations leads to a well-converged result. This may be 
compared with Fig. 2 of Cattaneo & Hughes ( 1996 ; hence- 
forth referred to as CH96), which suggests that, for a weak 
initital magnetic field, convergence is achieved only af- 
ter more than about 50 turnover times. However, there 
are important differences compared with their calculation. 
First, the dynamics in CH96 are very different, because 
the flow results from external forcing rather than moder- 
ately turbulent convection, and because the Lorentz force 
dominates, whereas it is negligible in our case. We have 
smaller fluctuations in the horizontal averages, partly be- 
cause the length scale of convection is always at least a few 
times smaller than the box size, whereas they are similar 
in CH96. Also, Re is smaller in comparison with CH96, 
and therefore the level of fluctuations. Note that in or- 
der to compare Re in Table [j] with the Reynolds number, 
U/fkf = 100, of CH96, where kf is the wave number of 
the forcing, Re should be multiplied with a factor of about 
1/27T. Secondly, the averaging in CH96 starts from an ini- 
tial state where there is no alpha effect, which exists only 
after 30 turnover times. This introduces a transient in the 
cumulative time average, which is absent in our case. As a 
result of these differences, we obtain convergence already 
if N « 6, even though this corresponds to a total averag- 
ing time, N ■ (t 2 — ii), of a few turnover times only. Fig. ||, 
top, shows mean time slices for the z-component of the 
EMF at different depths in the box. After an initial phase 
of approximately linear growth which, depending on the 
depth, lasts about 4 — 8 units, the EMF levels off to the 
equilibrium value. After that, it continues to evolve more 
slowly as it adapts to the changing mean magnetic field in 
accordance with Eq. The lower figure shows B x as a 
function of depth at consecutive moments. Starting from 
the initial condition, the mean magnetic field evolves while 
developing increasing gradients, especially near the upper 
and lower boundaries of the convective layer. In the time 
interval [5, 10] the magnitude of the gradients in the bulk 
of the unstable layer is d/dz\nB x w 0.5, and in the inter- 
val [20, 25] it is about 1. It follows that the best available 
time interval for calculating dynamo coefficients from run 
A45X is roughly between 5 and 10. 

As soon as the EMF levels off to its equilibrium value 
(||), the mean magnetic field already begins to develop 
considerable gradients at some depths. Nevertheless, the 
results are less contaminated with higher-order contribu- 



tions than those obtained by performing a simple time 
average accross a long simulation in which the magnetic 
field was initialized only once. An indication of the error 
in the coefficients is provided by the variation of the coef- 
ficients within the optimal time interval as shown by the 
shading in Figs. [| || In order to assess the possible dif- 
ference between the results of the improved method and 
those of the less accurate method of Paper I, we have 
in a few cases also calculated using the time interval 
[h,t 2 ] = [20, 25]. Although this choice of [h,t 2 ] yields the 
correct sign and order of magnitude, the higher-order con- 
tributions can amount to 40% in some cases. All results 
presented below were therefore derived using the earliest 
possible interval, so that the influence of the higher-order 
contributions is minimized. As is shown in Fig. |3j, bottom, 
the gradients of the mean field can be particularly strong 
near the top and the bottom of the unstable layer, so that 
at these locations the results should be treated with some 
caution. To avoid misunderstandings, we emphasize that 
the accuracy discussed here is distinct from the error in- 
troduced if the time sequence is too short to average out 
the fluctuations. 

4. Results 

The numerical simulations are grouped together in sev- 
eral series named A-E, each of which may contain several 
experiments performed at different latitudes and having 
different orientations of the imposed field. Each run is re- 
ferred to by a code such as A45X, which identifies it as 
belonging to series A, having 9 = 45° and an imposed 
magnetic field in the x-direction. The parameters of each 
series are summarized in Table [l]. They may be charac- 
terized as follows: series A is moderately turbulent and 
has a rotation rate appropriate for the base of the solar 
convection zone. Series B has the same parameters as A 
except for rotation, which is slower. Series C has the high- 
est Rayleigh number and the most turbulent convection. 
Series D has a Rayleigh number just above marginal sta- 
bility. Accordingly, the flow is laminar, and the Reynolds 
number is small. It has the slowest rotation of all. Series E 
is identical to D, except for the rotation rate. In series A, 
9 was varied between 0° and 90° with a stepsize of 15°; all 
other simulations were done for 9 = 45°. The optimal time 
interval, [tijfe], was determined to be [5, 10] for series A 
and B, [4,7] for C, and [15,20] for D and E. The value 
of N was in the range between 6 and 21. Although the 
total amount of time used in the averaging, N ■ (t 2 —tx), 
corresponds to several turnover times only, the combined 
spatial and temporal averaging as defined by Eq. |J) was 
in most cases sufficient to yield significant results. This is 
illustrated by the convergence of the cumulative average 
shown in Fig. ^, by the well-defined shapes of the curves 
of £ z and B x in Fig. |^, and by the narrowness of most of 
the shaded areas in Figs. 

In the following subsections we present the main re- 
sults of the simulations. As will become apparent, many 
of the results are new in the sense that apart from the 
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Fig. 4. Pumping velocity 7 (Eq. |4|) measured in units of 
0.01^/dg as a function of depth in units of d for series A. 
Each row, corresponding to a different value of 9, shows 
the three Cartesian components of 7 for the run indicated 
on the left. For the parameters, see Table 0. The black 
curves are time averages over the optimal interval as de- 
fined in Section 3. Shaded areas indicate the range of val- 
ues attained within the optimal time interval. The vertical 
scale is the same in each row, but differs from column to 
column. In series AOO only 7 Z was determined, while the 
other components of 7 are known to vanish on symmetry 
grounds. On the same grounds one knows that 7^ = in 
series A90. This was not imposed, so that this coefficient 
gives an indication of the error of the method. 



pumping effect in the vertical direction, none of the ef- 
fects have been numerically investigated before as far as 
we know. Also, this is the first numerical investigation in 
which pumping effects are quantified in terms of gamma 
coefficients. For that reason, and in order to compare the 
various pumping effects, it is worthwile to consider also 
the vertical pumping effect, although it has already re- 
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Fig. 5. Pumping velocity 7 (Eq. ||) measured in units of 
0.01^/dg as a function of depth in units of d from four 
different runs performed at 9 = 45° (cf. also Fig. || and 
Table 13). 



ceived much recent attention in the literature. Results for 



the coefficients a%j of all series are shown in Figs. |4|-g| in 
the form of the three Cartesian components of 7 and the 
independent elements of oq^ . From these, all components 
of ay can, if necessary, be derived using Eqs (§) and (§). 
The black curves in Figs. |4] |] are averages of a given co- 
efficient over the time interval [tijtg]- The shaded areas 
indicate the range of values attained at individual times 
within the interval \t\, £2], and are meant to give an error 
indication. Some confirmation that the error estimates are 
consistent, provided that they are used conservatively, can 
be seen in series A90, where 7 X , aty Z , a xx , a yy and a zz are 
known to vanish on symmetry grounds (Figs. ||, [3] and 
0). For these coefficients the shaded area is not always 
consistent with zero. This departure gives an additional 
estimate of the expected error. We reiterate that results 
for aij near the top and the bottom of the unstable layer 
should be treated with caution because of the enhanced 
possibility of contamination with higher-order terms. 



4.1. Latitudinal pumping effect 

In the latitudinal direction, the simulations yield a weak 
net general advection, measured by 7 X , which is equator- 
ward (73; > 0) in the bulk of the unstable layer and can 
be poleward near the upper surface and in the overshoot 
layer (Figs. |]-§|). At the pole and at the equator the effect 
vanishes; it is strongest near 9 = 60°. The symmetric con- 
tribution a yz is about as weak as 7^, but shows little sys- 
tematic behaviour as a function of 9 (Fig. ||). Both increase 
with increasing rotation rate, since they are rotationally 
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Fig. 6. Off-diagonal symmetric components of the ex- 
tensor (Eq. |^) measured in units of 0.0l\/dg as a func- 
tion of depth in units of d (cf. also Fig. and Table |l|). 
On symmetry grounds one knows that <J~ Z = in series 
A90. This was not imposed, so that this coefficient gives 
an indication of the error of the method. In series B45 and 
D45, which are characterised by rather slow rotation, a~. y 
is also not significant. 



induced (compare D45 and E45). As can be inferred from 



the relations 7^ = j x + a yz and 



lx 



yz 



|6j), the y and z components of the mean field have differ- 
ent latitudinal velocities (Fig. ||). In all simulations, B v 
tends to be pumped towards the equator (7^ > 0) in 
the bulk of the unstable layer. The situation with regard 
to the component B z is less clear. It is generally pumped 
towards the pole (7^ < 0) in the lower part of the un- 
stable layer, but the behaviour in the bulk of the unstable 
layer appears to vary with latitude. The average over the 
unstable layer turns out to be mostly positive (equator- 
ward) for both components (Fig. 0). The lowest-order con- 



tribution to latitudinal pumping in the expansion of the 
EMF presented by, e.g., Krause & Radler ( [L980D , p. 193, 
can be written as Aa^ = Afliejkigk^i + BQj€ikigk^U 
where A and B are true scalars, i.e. even functions of 
cos (9. Thus, the latitudinal pumping is at least of sec- 
ond order in Q. One mechanism that may explain the 
effect is provided by the rotational anisotropies studied 



by Kichatinov (1991). Indeed, the general tendency of his 
result is that the toroidal magnetic field, which for an ax- 
isymmetric mean field corresponds to our B y , is pumped 
towards the equator, while the poloidal field (i.e. B x and 
B z ) is pumped towards the poles (cf. his Eq. 3.26). 

4.2. Longitudinal pumping effect 

Surprisingly, the longitudinal pumping effect had so far 
not been identified in simulations of magnetoconvection, 
although it is often the strongest pumping effect except if 
the box is situated near the pole (Figs, column 2). 
The predominantly negative sign of j y signifies that the 
mean field is advected in the retrograde direction within 
the bulk of the unstable layer and the overshoot layer. 
In most cases, there is a thin layer near the top of the 
box where the field is pumped in the prograde direction. 
Unlike the vertical pumping effect (see below), here the 
symmetric part of ex-tensor dominates. From the relations 



(x) 

Tv = 



1v 



and 7y 



(*) _ 



7 y + a s xz (Eq. |)_ it follows 
that the longitudinal pumping velocities of B x and B z 
can be opposite. Since a xz is positive in the bulk of the 
unstable layer (Fig. ||, column 3), one finds that B x is 
pumped mainly in the retrograde direction, but B z mainly 
progradely near the top of the unstable layer, while rctro- 
gradely near the bottom. This behaviour is reflected in the 
volume averages (Fig. |). The longitudinal pumping effect 
is strongly dependent on latitude; it vanishes at the pole 
and peaks at the equator. The lowest-order term which 
contributes to longitudinal pumping (Krause & Radler 
1980, p. 193) is Actij = Ag i ^lj+Bgj^l il where A and B are 



true scalars. Clearly, the effect requires both gravity and 
rotation. The antisymmetric part, corresponding to gen- 
eral longitudinal pumping, has been computed by Riidiger 
& Kichatinov (1993) in the form of a^(Q,iGj — iljGi), 



where the coefficient proved to be positive-definite. 
Here G is the gradient vector of the turbulence intensity, 
i.e. G = Vlnw rms . One obtains 



(Eq. 7y oc rising- 



,du 2 
dz 



while for the vertical pumping the simple relation 



7z °c 



du 2 
dz 



(14) 



(15) 



results, with only a very weak dependence on the Coriolis 
number. Thus j y and ^ z should always satisfy the relation 
7y7z < 0, which is indeed the case in all the presented 
simulations. According to Eqs (14) and (|l5|), the nodes 
which exist in all profiles of 7 a and 7 2 should be very 
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Fig. 7. Diagonal components of the a-tensor (Eq. mea- 
sured in units of 0.01 \fdg as a function of depth in units 
of d from runs of scries A performed at different values of 
9 (cf. also Fig. ^ and Table 0). On symmetry grounds one 



knows that a 2 



in series A90. This was 



not imposed, so that these coefficients give an indication 
of the error of the method. 



close together, which is really true in Figs. |4|-|3[ For com- 
pleteness we also note that longitudinal pumping cannot 
be attributed to the t ype of rotational anisotropics studied 
by Kichatinov ( 1991 ), since the resulting pumping effects 
exist only in the radial and latitudinal directions. 

4.3. Vertical pumping effect 

The main features of the vertical pumping effect are as 
follows. The direction of vertical pumping is downward 
(72 > 0) in the bulk of the unstable layer and in the over- 
shoot layer below it (Figs. Near the top of the box 
there is a thin layer where the pumping is directed up- 
wards (-f z < 0). The strongest effect is observed in the 
laminar cases (D45 and E45), where upward and down- 
ward speeds can amount to about 15% of the rms convec- 
tive velocity (Table |l|). In the more turbulent runs, the 
maximum speeds are smaller. Since the top layer where 
the pumping is upward is less pronounced in the turbulent 
runs, the decrease is less apparent in the volume-averaged 
coefficient (Fig. ^) because there is less cancellation. There 
is little dependence on rotation (compare A45 with B45 
and D45 with E45) or latitude (Figs. || and ||). Also the rms 
convective velocity has a depth dependence that is similar 
in all runs, namely a plateau at the top of the domain 
and below that a monotonic decrease with depth. This in 
combination with the robustness of the vertical pumping, 
in particular its virtual independence of rotation, suggests 
an explanation in terms of 'turbulent diamagnetism', i.e. 




Fig. 8. Diagonal components of the a-tensor (Eq. ||) mea- 
sured in units of 0.01^/dg as a function of depth in units 
of d from four different runs performed at 9 = 45° (cf. also 
Fig. I and Table @). 



by transport down the gradient of the rms convective ve- 
locity (Eq. |l5|) . Analytical expressions for the diamagnetic 
pumping effect derived under the FOSA assumption were 
given by Kichatinov & Riidiger (1992). 

If the Rayleigh number is increased from 10 6 (series B) 
to 4 • 10 6 (series C), so that the convection becomes more 
strongly turbulent, \ j z | decreases. This would appear to be 
at odds with the result of Tobias et al. ( |2001| ) , who found 
more efficient pumping for increased Rayleigh numbers. 
But they measure the pumping effect through the ratio of 
magnetic fluxes in the unstable layer and the overshoot 
layer. This ratio is determined both by transport proper- 
ties and storage properties of the magnetic field, i.e. by 
the balance between pumping, turbulent diffusion in the 
convection zone, and molecular diffusion in the overshoot 
layer (ij). If 77 is decreased, flux concentrations decay more 
slowly, so that the storage capacity is increased. Since the 
increase of the Rayleigh number in Tobias et al. ( [2001 ) is 
accompanied by a decrease of 77, the flux ratio increases, 
even if the pumping is somewhat less efficient. 

In series A it can be seen that the symmetric contribu- 
tion a^y is positive and that it increases if 9 increases from 

-at% and 



15° to 75° (Fig. §). Using the relations 7, 



(x) 



a xy (Eq. ^) one can infer from this that the 



(y) 

pumping speeds for B x and B y deviate increasingly from 
one another: the downward pumping of B x becomes faster, 
and the top layer where B y is pumped upward becomes 
deeper; this is reflected in Fig. ||. A comparison of runs 
B45 with A45 and D45 with E45 shows that ecL increases 

y 

with increasing rotation rate, while j z is hardly affected 
(Figs. In run D45, which has the weakest rotation of 
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Fig. 9. Averages over the unstable layer (0 < z < 1) of 
mean-field transport coefficients for series A (Figs. 
measured in units of 0.01^/dg, as a function of angular 
distance to the south pole. The field-direction dependent 
pumping effects, shown for clarity in the bottom six fig- 
ures, are related to the other coefficients through Eq. @. 
A rough indication of the errors is provided by the shaded 
areas in Figs. |]-|]; the resulting error in the z-averages is 
typically 0.15 in units of O.Olv^- Thus the results for a yz 

are not significant, while those for "f x , a zz , and ji^ 
marginally significant. 



arc 



all, the value of a xy is not above the noise level, and the 
same is probably the case for run B45. From the symmetry 



considerations of Krause & Radler (1980), p. 193, it is in- 
ferred that the lowest-order contribution to j z stems from 



the term Aa.-, = Ae; 



while that to a™ stems from 



Actij = AQiEjkigkfli + BQje lk ig k Qi, where A and B are 
true scalars. The latter term also gives rise to latitudinal 
pumping, as mentioned above. Unlike j z , the coefficient 



a xy is rotationally induced, and this is confirmed by the 
simulations. 

4.4. Alpha effect 

As a byproduct of the determination of the pumping ef- 
fects, we obtain results for the diagonal part of the a- 
tensor. These results are an extension to what was re- 
ported in Paper I, where all simulations were done for 
Ra = 50 000 and 9 = 0°. The sign and magnitude of ot yy 
are of special interest, since this component is responsi- 
ble for the generation of the large-scale magnetic field in 
solar-type stars. The depth dependence of a xx and a yy 
has a common shape in all series, provided there is suffi- 
cient rotation, namely a negative sign in the bulk of the 
convection zone, and a positive sign in the overshooting 
layer (Figs. These features, which are expected for 

the southern hemisphere, are explained, e.g., in Paper I. 
Series A (Fig. [?]) provides detailed results on the latitude 
dependence. The amplitudes of a xx and a yy increase with 
angular distance from the equator up to a point close to 
the south pole, more or less consistent with the commonly 
assumed cos ^-function (Fig. |J). At the pole, there appears 
to be an unexpected dip. The amplitudes of a xx and a yy 
are approximately the same near the pole - for 8 = 0° they 
are equal on symmetry grounds - but closer to the equator 
a yy is somewhat larger (especially E45; Fig. ||). Without 
rotation, the diagonal a-terms vanish (the results in Fig. |^ 
for D45, which has the weakest rotation, are probably not 
significantly different from zero), and the increase of a xx 
and a yy with rotation (compare B45 with A45 and D45 
with E45 in Fig. ||) confirms the results of Paper I. The 
same can be said for a zz , which for Co ^ 2 has the op- 
posite sign of a xx (tentatively visible for B45 in Fig. || 
but clear in Paper I) . If rotation increases beyond a cer- 
tain point, a zz as a function of depth developes multiple 
sign changes (cf. series A), and its amplitude falls behind 
that of a xx , even more so if one also averages over depth 
(Fig. ^). This is the rotational quenching of the vertical 
alpha effect reported in Paper I. 



5. Discussion 

We have performed numerical simulations of magnetocon- 
vection in a Cartesian box with a rotating frame of refer- 
ence which is meant to represent a section from the lower 
part of a stellar convection zone. Through an inversion 
of the relation between the mean magnetic field and the 
electromotive force, we have measured pumping effects in 
all three directions as well as the symmetric elements of 
the alpha tensor. For the first time, the existence of non- 
radial pumping effects is hereby established numerically, 
and quantified in terms of gamma coefficients. Previously, 
their existence was known only on the basis of symmetry 
considerations in the fram ework of mean-field electrody- 
namics (Krause & Radler 1980 , p. 193), and analytical 



results base d on t he first-order smoothi ng app roximation 
(Kichatinov |l991 ; Riidiger & Kichatinov 1993 ). Such ana- 
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lytical results are available only for some of the coefficients 
that appear in the expansion of the EMF given by Krause 
& Radler ( 1980). Moreover, the validity of results based on 
FOSA is questionable for real and simulated stellar con- 
vection zones alike, and they can be readily evaluated only 
for simple turbulence models but not for actual convective 
flows. For these reasons, it is difficult at present to explain 
all of the phenomena that are found in the simulations, al- 
though we did observe that several qualitative features of 
the expressions derived using FOSA are confirmed by the 
simulations, as was the case in Paper I. 

The Rayleigh number was varied between 10 , corre- 
sponding to laminar convection, and 4 • 10 6 , correspond- 
ing to mildly turbulent convection. Within this parameter 
range, we found remarkably consistent behaviour of the 
dynamo coefficients. Clearly, the regime of stellar convec- 
tion, which is characterized by Rayleigh numbers many 
orders of magnitude larger, remains far out of reach be- 
cause of the limited spatial resolution. The main resulting 
difference between the simulations and stellar convection 
is as follows. In the lower part of stellar convection zones, 
the superadiabaticity is very small (typically 5 « 10 -6 ) 
and, according to the standard mixing-length expressions 
(e.g., Stix 1989, p. 198), this results in velocities of the 
1 or Ma ss 10~ 4 . In the simulations, con- 



up to about 0.1u rms in the turbulent runs for 9 = 45° 



while the volume averages (7: 



(y) 



a 



in Fig. 



are 



Stix 
order of 50 m s 



vection is less turbulent. After the onset of convection, the 
superadiabaticity in the unstable layer decreases from its 
initial value, 82 — 0.1, to a value of at least about 10 -2 
(series C). The convective velocities are therefore much 
larger than in the Sun, if expressed in m s _1 or as a Mach 
number. We may illustrate this for series A, where the 
unstable layer covers 3.7 pressure scale heights (Table |l|). 
Counting upwards from the bottom of the solar convection 
zone, this corresponds to about 90 000 km. Setting g « 540 
m s~ 2 , one would obtain \fdg~ « 2 • 10 5 m s -1 for the unit 
of velocity. Thus the rms convective velocity of series A 
(Table |l|) would correspond to about 2 • 10 4 m s^ 1 , i.e. ap- 
proximately 400 times larger than the solar value. About 
the same ratio holds for the Mach numbers. Only simu- 
lations with a much higher resolution and much smaller 
entropy gradients can improve this ratio. However, the ro- 
tation rate was chosen such that the Coriolis number, i.e. 
the relative importance of the Coriolis force in the equa- 
tion of motion, is within the appropriate range. Also, it 
is to be expected that beyond some Rayleigh number, a 
further decrease of the molecular viscosity or the magnetic 
diffusivity affects only the smallest scales while having lit- 
tle or no effect on the large scale. From the consistency 
of the results one may therefore suspect that the current 
simulations already provide some clues about the mean- 
field transport properties of stellar convection, if properly 
scaled. 

In order to assess the astrophysical relevance of the re- 
sults, the pumping speeds should therefore be compared 
with the rms convective velocity for the unstable layer 
(Table [l]) rather than be directly expressed in m s _1 . 
The latitudinal pumping effect for the toroidal mean field 
(By) was found to be mainly equatorward with speeds 



less than about 0.03u rms . The equatorward motion in the 
solar butterfly diagram corresponds to a speed of about 
1 m s _1 , so that one may conclude that the latitudinal 
pumping effect of the toroidal mean field has the same or- 
der of magnitude. We therefore point out the possibility 
that the equatorward motion of solar magnetic activity 
is not the result of a dynamo wave or meridional bulk 
motion, but of latitudinal pumping of the toroidal mean 
magnetic field. Results for the latitudinal pumping of the 
radial mean field do not show consistent behaviour; the 
poleward motion at the top of the box found in some runs 
amounts to at most 0.1u rms , but this should be treated 
with caution because of contamination with higher-order 
terms due to the strong gradients in the mean field near 
the top of the box. Since the simulations are not meant 
to represent the photosphere, the relevance of this result 
for the observed poleward motion of poloidal fields at the 
solar surface is in any case hard to establish. We merely 
point out the possibility that the observed motion is the 
result of pumping of the deep-seated global poloidal field 
rather than of meridional bulk surface motion. The longi- 
tudinal pumping attains maximal speeds of about 0.2u rms 
in the retrograde direction. Again, the mainly prograde 
pumping of the poloidal mean field near the top of the 
box should be treated with caution because of the pos- 
sibility of strong higher-order contributions. Perhaps, if 
carried over to the surface of the Sun, it could explain 
the observed motion of the active longitudes in the so- 
lar wind relative to solar surface rotation (Ruzmaikin et 
al. 2001). We also note that preliminary calculations of a 
mean-field model indicate that the inclusion of longitudi- 
nal pumping in an a 2 -type dynamo can lead to oscillating 
solutions, which are rare for such models. This can be rel- 
evant for explaining magnetic cycles in stars which do not 
have differential rotation (Riidiger et al. 2002 ) . The down- 
ward pumping in the turbulent runs attains speeds of up 
to about 0.1u rms , and the volume average about half of 
that. In the lower part of the solar convection zone this 
translates to a few m s . Such speeds may be sufficient to 
overcome the magnetic buoyancy of flux generated in the 
convection zone and transport the flux downward into the 
overshoot layer. This layer coincides more or less with the 
tachocline, a region where the rotation rate changes rather 
suddenly with depth. Through stretching of field lines, the 
shearing motions can generate from poloidal field compo- 
nents a strong field of predominantly toroidal orientation, 
as is required in order to account for various features of 
the large-scale solar magnetic field. 

The origin of the poloidal magnetic field and its locus 
of generation still remain an open question. Accordi ng to 



one likely scenario, th e inte rface dynamo (Parker |1993 



Ossendrijver & Hoyng 1997, Charbonneau & MacGregor 



1997 ), the poloidal field is generated by an alpha effect in 
the lower half of the convection zone. Our results for the 
coefficient a vv , which corresponds to in spherical coor- 
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dinates, are of particular interest for interface-type stellar 
dynamo models. We should stress however that it cannot 
be decided from the present simulations whether the large- 
scale magnetic field of solar-type stars is generated by the 
convective alpha effect investigated here or by a different 
alpha effect such as that resulting from the buoyancy in- 
stability of magnetic flux tubes in the overshoot layer. We 
confirm that the convective alpha effect is strongest near 
the poles. For the Sun, this should, in combination with 
the differential rotation near the poles, as inferred from he- 
lioseismic inversions, result in strong dyna mo action near 
the poles (Rudiger & Brandenburg 1995 ), which is not 
observed - sunspots are concentrated in belts that extend 
to a latitude of at most 40° on either side of the equa- 
tor. This well-known problem may perhaps be resolved 
by taking into account meridional circulation, and by not- 
ing that the buoyancy instability of toroidal magnetic flux 
tubes anchored in the overshoot layer is reduced or sup- 
pressed at high latitudes (Ferriz-Mas & Schussler 1995| ). 
The sign of a yy confirms what is generally found in MHD 
simulations (e.g. Paper I), namely a negative sign in the 
bulk of the convection zone on the southern hemisphere, 
i.e. a positive sign on the northern hemisphere, and a sign 
change at the bottom of the convection zone. This dis- 
tribution of signs puts a tough constraint on theories of 
the solar dynamo. Generally, mean-field dynamo models 
that incorporate solar differential rotation, in particular 
a positive dtt/dr near the equator near the base of the 
convection zone, produce a butterfly diagram with the 
required equatorward migration of the magnetic activ- 
ity belts only in combination with a negative on the 
northern hemisphere. There are various ways in which this 
dilemma could be resolved; perhaps the relevant dynamo 
effect is located at the very bottom of the convection zone 
where the sign has changed. If is positive on the north- 
ern hemisphere, equatorward migration can be obtained 
if a suitable meridional circulation is included (Choudhuri 



et al. |1995j Durney [1996; Kuk er et al. gOOlf Rudiger et 
al. |200l|; Rudiger & Arlt 2002|). As mentioned above, the 



present results suggest the alternative possibility of lat- 
itudinal pumping. This suggestion needs however to be 
verified in mean-field models that take into account all 
relevant dynamo coefficients as well as solar differential 
rotation. 
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